race = unique(analysis_df$race_recode)
foreign = unique(analysis_df$birthplace_recode_first)
occ = unique(analysis_df$occ_recode)


first_stage_results = NULL
for(i in race){
  for(j in foreign){
    for(g in occ){
      
      lm1 <- try(lm_robust(as.formula(paste0('govt_', i, '_x_', j, '~ treat + ', i, '_x_', j, '+ govt_', i, '_x_', j, '_lag + ', i, '_x_', j, '_lag')),
                           data=analysis_df_3 %>% 
                             filter(occ == g) %>%
                             group_by(city) %>%
                             arrange(city, YEAR) %>%
                             mutate_at(vars(govt_occ_total:white_x_native_born), list(lag=lag)), 
                           se_type='stata',
                           clusters = city), TRUE)
      
      if(class(lm1) == "try-error"){next}
      
      first_stage_results <- first_stage_results %>%
        bind_rows(unlist(c(tidy(lm1)[2,], race = i, foreign = j, occ = g)))
    }
  }
}

top_results <- first_stage_results %>%
  as.data.frame() %>%
  mutate_at(vars(estimate:df), function(x){as.numeric(as.character(x))}) %>%
  filter(outcome %in% c("govt_black_x_native_born", "govt_white_x_foreign_born", "govt_white_x_native_born")) %>%
  mutate(plot_x = case_when(race == "black" ~ "black",
                            race == "white" & foreign == "foreign_born" ~ "foreign born white",
                            race == "white" & foreign == "native_born" ~ "native born white")) %>%
  mutate(occ = str_replace(occ, "\\_", " ")) %>%
  select(estimate, std.error, occ, plot_x, p.value)




foreign_orig = unique(analysis_df$birthplace_recode_foreign)
occ = unique(analysis_df$occ_recode)

first_stage_results = NULL
for(i in foreign_orig){
  for(g in occ){

    lm1 <- try(lm_robust(as.formula(paste0('govt_', i, '~ treat + ', i,'+ govt_', i, '_lag + ', i, '_lag')),
                         data=analysis_df_4 %>% 
                           filter(occ == g) %>%
                           group_by(city) %>%
                           arrange(city, YEAR) %>%
                           mutate_at(vars(govt_occ_total:russian), list(lag=lag)), 
                         se_type='stata',
                         clusters = city), TRUE)
    
    if(class(lm1) == "try-error"){next}
    
    first_stage_results <- first_stage_results %>%
      bind_rows(unlist(c(tidy(lm1)[2,], race = i, occ = g)))
  }
}

  
bottom_results =
  first_stage_results %>%
  as.data.frame() %>%
  filter(race != "other") %>%
  mutate_at(vars(estimate:df), function(x){as.numeric(as.character(x))}) %>%
  mutate(plot_x = race) %>%
  mutate(occ = str_replace(occ, "\\_", " ")) %>%
  select(estimate, std.error, occ, plot_x, p.value)
  

results =
  bind_rows(top_results, bottom_results) %>%
  arrange(desc(plot_x))


table <-
  results %>%
  filter(occ %in% c('blue collar', 'white collar') & std.error!=0) %>%
  mutate(occ = str_to_title(occ),
         group = str_to_title(plot_x),
         group = case_when(group=='Foreign Born White' ~ 'Foreign W.',
                           group=='Native Born White' ~ 'Native W.',
                           group=='Black' ~ 'Native B.',
                           T ~ group)) %>%
  mutate(p.value = round(p.value, 3)) %>%
  select(Group = group, Occupation = occ, Estimate = estimate, SE = std.error, `P-Value` = p.value)

table <- stargazer(table[16:1,], title = "Effect of Civil Service Reform, Lagged-DV Models", summary=F, rownames=F)

#cat(table, sep = '\n', file='../Apps/Overleaf/merit paper/outputs/lag_dv.tex')
cat(table, sep = '\n', file='./replication_file/_5_outputs/tables/table_a4.tex')


